Giant Charge Inversion of a Macroion Due to Multivalent 
Counterions and Monovalent Coions: Molecular Dynamics Study 
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We report molecular dynamics simulation of the (overall neutral) system consisting of an immobile macroion sur- 
rounded by the electrolyte of multivalent counterions and monovalent coions. As expected theoretically, counterions 
adsorb on the macroion surface in the amount much exceeding neutralization requirement, thus effectively invert- 
ing the sign of the macroion charge. We find two conditions necessary for charge inversion, namely, counterions 
must be multivalently charged and Coulomb interactions must be strong enough compared to thermal energy. On 
the other hand, colon condensation on the multivalent counterions similar to Bjerrum pairing is the major factor 
restricting the amount of charge inversion. Depending on parameters, we observe inverted charge up to about 
200% the original charge of the macroion in absolute value. The inverted charge scales as ~ C}^'^ when C < 1 3'iid 
crosses over to ~ for C, > 1, where C, = (Ao/rs)^, is the Debye screening length in the electrolyte and Aq is 
the distance between adsorbed counterions under neutralizing conditions. These findings are consistent with the 
theory of "giant charge inversion" [Phys. Rev. Lett., 85, 1568 (2000)]. 

PACS numbers: 61.25.Hq, 61.20.Qg, 87.14.Gg, 82.70Dd 



I. INTRODUCTION 

Correlation effects in the systems of charged particles, 
such as plasma or electrolyte solution, are well known 
since the works by Debye and Hiickcl in 1923 0. Clas- 
sical intuition suggests that correlation can be viewed as 
screening in which a cloud of ions around, say, positive 
particle is slightly dominated by negative counterions, 
such that for an outside observer (who measures the elec- 
tric field) the shield of predominantly negative charges 
effectively reduces the central positive charge. Recently, 
a significant attention has been attracted by the notion 
that much more dramatic effect is possible in the system 
with strongly charged ions [?] . Namely, instead of charge 
reduction due to the shielding, it is possible to observe 
charge inversion due to the " over-screening" . Further- 
more, it was shown a year ago that the inverted charge 
may be quite large, even larger in absolute value than the 
original bare charge, giving rise to the concept of "giant" 
charge inversion |^ . 

In the present paper, we use molecular dynamics simu- 
lation technique to address the question of possible limits 
of charge inversion. Overall, we confirm the theoretical 
prediction Q and observe " giant" charge inversion, with 
the ratio of inverted and bare charges reaching up to 
about 1.6 (in absolute value). 

Although we consider here only primitive schematic 
model with spherical ions immersed in the medium of a 
constant dielectric permeability e, this should be viewed 
as the step towards better understanding of such first 
magnitude scientific problems as, e.g., that of chromatin 
structure. Indeed, chromatin represents a complex of 
strongly negatively charged thread of DNA with posi- 
tively charged smaller protein molecules. For instance, 
virtually every paper on charge inversion mentions the 



fact that protein core of a nucleosome particle ||lj| car- 
ries lesser amount of positive charge than the amount of 
negative charge on the wrapped around DNA. On a sim- 
pler level, complexes of polycations and polyanyons were 
under scrutiny for a long time fl^ ], as well as complexes 
of charged polymers with charged colloids [p^ . 

In theoretical aspect, the most advanced treatment of 
charge inversion is due to Shklovskii and his co-workers 
In these works, the universal physical mechanism 
behind charge inversion is recognized as correlations be- 
tween shielding ions (see also brief review article Q ) . In 
the work |^ , the idealized image of these shielding coun- 
terions forming a Wigner crystal on the surface of the 
shielded macroion was emphasized (see also earlier work 
[Q). However, it was mentioned in and addressed in 
more details in that in most real cases, correlations 
are not quite as strong as to produce a crystal, but suffi- 
cient to maintain short range order, and, therefore, cor- 
relation energy is similar to that of a crystal. Obviously, 
this mechanism is operational when shielding ions are 
strongly charged. Furthermore, it was realized that the 
best situation for charge inversion occurs when monova- 
lent salt is present in addition to strongly charged ions j|] 
(see also [^). Salt ions, as their charges are small, behave 
in a " traditional" way; they simply screen all interactions 
at the distance about Debye length Vg. However trivial 
itself, this leads to a dramatic increase of charge inver- 
sion, because the attraction of a counterion to its Wigner- 
Seitz cell on the macroion surface is over a significantly 
shorter range than the repulsion of a counterion from 
the uncompensated charge of all other counterions. For 
completeness, we mention here also recent works devel- 
oping charge inversion theory to include polyelectrolyte 
ions [pl pdl , as well as more formal theoretical approaches 

Charge inversion has been seen several times in simu- 



lations, starting from the pioneering work |17| . In recent 
works Jl8|- p3| , p^ computer simulations were reported 
along with various ways to re-derive and re-examine the 
concept of lateral correlations between counterions as the 
driving force behind charge inversion. The authors of 
| p2| , pl| reported quite impressive agreement between the- 
oretical conjectures and their computation data. How- 
ever sophisticated, these simulations concentrated on the 
cases of no added salt and of abundance of counterions. 
In other words, they only examined the very dilute ex- 
treme with respect to macroions assuming at the same 
time finite concentration of counterions. Our first intent 
in the present work is to relax this serious restriction 
and to simulate a realistic model in which thermody- 
namic cost of adsorption of counterions on the surface 
of a macroion is contributed by both the events on the 
macroion surface and in the surrounding solution. 

The other closely connected goal of our present study 
has to do with the following delicate aspect of the " giant" 
charge inversion scenario. In order to make correlations 
and charge inversion stronger, one is tempted to choose 
larger ratio of Coulomb energy to thermal energy. But 
when it is too large, the small salt ions start to condense 
on the surfaces of counterions, effectively reducing their 
charge. Therefore, charge inversion is expected to be the 
strongest in the intermediate regime, when correlations 
between counterions are already strong, but condensation 
of small ions on them is still weak. Therefore, we want 
to check in the present work computationally how robust 
is this theoretical prediction. 

To achieve the above stated goals, we perform in the 
present paper molecular dynamics study of the system 
consisting of a single macroion, large number of multiva- 
lent counterions, and a multitude of monovalent colons 
immersed in a Langevin fluid. It is worth noting that 
hydrodynamic effects, which may be of significant im- 
portance for interactions between colloidal particles away 
from thermodynamic equilibrium p^-p^, are totally ig- 
nored in the present study, because we concentrate on 
the equilibrium aspects only. 

The paper is organized as follows. The simulation 
method and parameters are described in Section |l[ In 
Section III, by direct measurement of the peak height of 
the radial charge distribution we show that giant charge 
inversion takes place when the following two conditions 
are simultaneously met: (1) multivalent counterions with 
valence Z >2 are present, and (2) Coulomb energy pre- 
vails over the thermal energy at the length scale of a 
single ion size, a: T — Z'^e^ / eaksT > 1. We study in 
details the dependence of charge inversion on the radius 
and charge of the macroion, the valence and density of 
counterions and colons, and temperature. For large den- 
sity and valence of counterions, the amount of inverted 
charge increases linearly with ionic strength, and reaches 
up to 200% the original macroion charge. 

Extension of the present work to the case under elec- 



trophoretic environments is discussed in a separate paper 
p7t in which the effect of an applied electric field on the 
charge inversion process is investigated with the use of 
molecular dynamics simulation. 



II. SIMULATION METHOD AND PARAMETERS 
A. Equations 

Specifically, we consider the following model. The sys- 
tem includes: a single macroion with negative charge 
Qq < 0, multivalent counterions with a positive 
charge Ze each, and N~ monovalent colons with a 
negative charge (— e) each (e > is the elementary 
charge). Overall charge neutrality is strictly enforced: 
Qq + ZeN'^ — eN~ — 0. All ions arc confined within 
the three-dimensional simulation domain having spheri- 
cal shape with radius Rm- The macroion is considered 
immobile; it is placed at the origin (center of the domain), 
and all other ions are mobile. All ions are supposed to be 
of spherical shapes, with macroion having radius Rq and 
all mobile ions having identical radius a; a serves also as 
a unit of length. 

The (classical) molecular dynamics simulation solves 
the Newton-Langevin equations of motion 



TO— = -V$(rO - V0(ri) - lya-Vi + Fth , 



(1) 



where the potentials and 0, describe interactions of a 
given ion with other mobile ions and with the macroion, 
respectively: 
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Here, and are the position and velocity vectors of the 
i-th particle, — jr^ — rj|, e the dielectric constant, eLj 
the Lennard-Jones energy. As regards boundaries, we as- 
sume elastic reflection every time when a mobile ion hits 
either the domain boundary at r — Rm or the macroion 
surface at r = Rq. The last two terms of Eq.(|l|) repre- 
sent the Langevin thermostat due to surrounding neutral 
medium. The Stokes formula for a sphere is adopted for 
the friction term with v being the friction constant, and 
Fth is the random (5-correlated thermal agitation. 

The inertia term is retained in the momentum equation 
for numerical stability of the electrostatic forces, masses 
of all mobile ions are assumed identical, equal to m. This 
leads to the choice oiuj^^ as the natural time unit, where 
ujp = (47rnoe^/eTO)^/^ is plasma frequency and uq the 
average ion number density. 



B. Parameters 

It must be born in mind that phenomena resembhng 
charge inversion may occur when other forces, apart from 
Coulomb electrostatic ones, operate in the system (in- 
cluding complicated helical shape of the molecules in- 
volved; see, for instance, [p8|). In this study we are in- 
terested in the situation when pure electrostatic forces 
dominate. Accordingly, we choose e^j = (l/12)e^/ea; 
this corresponds to the depth of Lennard- Jones potential 
well equal to — eLj/4 = — (l/48)e^/ea, which means that 
Lennard-Jones attraction force is very small compared to 
Coulomb force even at the distance of ion size a and even 
for monovalent ions. 

We also consider densities at which short range repul- 
sion (excluded volume effect) between ions is not impor- 
tant, as volume fraction of particles in the simulation 

\N+ +N-)/{Rl,-Rl), is 



domain, (p = 
small, about i 
this paper. 

By contrast. Coulomb interactions are strong. To be 
more specific, there are several relevant parameters con- 
trolling different manifestations of Coulomb forces. First 
of all, multivalent Z-ions are attracted to the macroion 
and can be adsorbed on its surface. This is controlled by 
the parameter 



in water (e « 80). In particular, for the typical small 
ions, for which a « 4A (counting attached water), we get 
Ta w 1.7ZA, which is roughly between 4 and 10 for Z 
between 2 and 7. As regards F = Fa (a/4i?o) \/ZQ/e, it 
may be greater than Fa if macroion is strongly charged 
(Q/e is large). 

Note also that under typical conditions, such as to « 
SOtojj and no ~ (1/loA)'^, where mu is proton mass, and 
ng the average density of counterions, the characteristic 
frequency and time are about ujp « 6.6 x lO^^s^^ and 
LU~^ « l.hps. 

In our molecular dynamics experiment, the initial po- 
sitions of CO- and counterions are distributed randomly 
between the two spheres Rq < r < Rm, each ions hav- 
ing the velocity that satisfies the Maxwell distribution. 
The integration of the equations of motion is done with 
the use of the leapfrog method which is equivalent to 
Verlet algorithm |29|. The time step of integration is 



0.05 or less for aU cases considered in = O.Olujp \ and simulation runs are executed up to 

SOOOw"^ at which time the peak height of the inverted 
charge Eq.(^) has become stationary. 

Below, in Section [11 , we report the simulation results 



To = 



ZeQ 
eRoT 



(3) 



(for the temperature T, we use energy units and omit 
Boltzmann constant fcs). Second, monovalent colons are 
attracted to the multivalent counterions (Z-ions), and 
can condense there, which is controlled by the parameter 



F„ = 



eaT 



(4) 



A little more delicate matter is the possible correlation 
between repelling ions, particularly those adsorbed on 
the macroion. This is characterized by F = Z'^e^/eAT, 
where A can be estimated as the distance between two 
adsorbed counterions in the situation when the num- 
ber of adsorbed counterions is just sufficient to neutral- 
ize the macroion, that is {Q/Ze)TT{A/2)'^ = AttRq or 
A = ARoy^Ze/Q. Thus, 



concentrating on the general properties of the charge in- 
version: its dependence on the radius and charge of a 
macroion, the valence and density of counterions, and 
temperature. While changing the parameters, the elec- 
trostatic binding energy of counterions to the macroion 
is kept constant by fixing Fq Eq.(^. 

In the present study, the following values of parame- 
ters are considered " standard" and used unless otherwise 
specified: radius of the macroion Rq = 3a, its charge 
Qo = — 28e (assumed negative), valence of the counteri- 
ons Z — 7, and the number of the counterions and colons 
= 52 and N~ = 336, respectively. The radius of the 
outer boundary sphere is Rm = 20a. The temperature 
is chosen such that Fq = 4.2Z. 

To support physical intuition, it is useful to estimate 
the Debye screening length. Naive application of stan- 
dard formula yields 



47re 



Z^N^ 



N- 



{RM/af - 



(Ro/af 



-1/2 



1/2 



3Fa(ZiV+ + Z-iiV-) 



(6) 



^3/2g3/2Ql/2 

AeRoT 



(5) 



In principle, there is also other similar F parameters 
which control correlations between various ions in the 
bulk; in this work we do not address this aspect. 

In the present study, we typically look at the Fq values 
in the range F = 6 ^ 80. For the estimates, it is useful 
to keep in mind that Bjerrum length £b = e^/eT is close 
to 7 Aunder typical conditions - at room temperature 



which is about 0.5a under the "standard" conditions. 
This result may seem surprising, as physically screening 
length cannot be smaller than the size of smallest ions 
pof . Of course, such a small value of screening length 
indicates very strong Coulomb interactions in the bulk 
solution. This fact can be also seen differently, by not- 
ing that the parameters controlling validity of the lin- 
earized Debye-Hiickel theory for the plasma away from 
macroion are Z(f)]^^Ta and (j)]/^Ta/Z , and they are both 
large compared to unity, about 10-200 (Z=3-7) and 1-4, 



respectively. (These parameters mean ratio of Coulomb 
between particles of respective signs and thermal energies 
at typical distances - controlled by densities.) Thus, we 
consider the conditions under which plasma outside the 
macroion is very nonlinear. Physically, this is manifested 
by extensive condensation of coins on multivalent counte- 
rions, as will be seen in the results below. Such condensa- 
tion is analogous to Bjerrum pairs formation. Condensa- 
tion means that effective charges of particles are reduced, 
and also the effective density of charged particles is low- 
ered. That leads to the increase of real screening radius 
which attains some respectable value. We do not attempt 
to estimate it, as we do not rely on any particular the- 
ory. Instead, we will just see what molecular dynamics 
show. We shall see that condensation of colons is the 
major factor limiting the extent of charge inversion. 

Another interesting quantity to estimate is Gouy- 
Chapman length associated with the surface of a 
macroion, 

2ZTiea rjQo/el ^' 

(where a — \Qo\ /AttR^) turns out to be about Q.lba/Z. 
(Strictly speaking, A is defined for the plane, not spherical 
surface; however, since \/Rq ~ 0.05 ^ 1, defining A 
based on plane geometry is reasonable.) 

For the standard run, it takes about 2.5 x lO'^w"^ be- 
fore a state is reached which can be assumed equilibrated, 
at least in terms of the inverted charge being stationary. 

III. SIMULATION RESULTS 
A. Observing charge inversion 

1. Standard regime 

The results of our simulations are presented in the 
Figures p]-p^. Figures ^ and || present typical results 
of runs performed under what we call "standard" con- 
ditions. Specifically, Fig. |^ shows a snapshot of the 
spatial distribution of counterions and colons around the 
macroion after charge distribution has become station- 
ary. Since our simulation includes hundreds of particles, 
it is impossible to "see" them in any meaningful way; 
what we can see, however, is the configuration of ions in 
the immediate vicinity of the macroion surface. This is 
shown in Fig. ^ in which only the ions residing in the 
thin layer Rq <r < Rt^ + ia are depicted. 

As seen in Fig. [l], counterions (light blue) attach right 
on the surface of the macroion with a lateral spacing, 
while colons (dark blue) stay some distance away from 
the macroion surface. It is clear that lateral correla- 
tions are present between counterions, particularly be- 
cause there are no pairs in which counterions are close to 



each other. Not surprisingly, however, this correlations 
are much weaker than in the case without colons exam- 
ined in although counterions are correlated in Fig. 
|l|, their spacings are not regular and cannot be identified 
as Wigner crystal. As regards colons, they are seen to 
condense on the top side of the counterions, presumably 
because of strong repulsion of the colons from macroion 
surface. We note here that this condensation of colons 
on the counterions is responsible for limiting the amount 
of charge inversion. In the configuration shown in Fig. 
|l|, the numbers of counterions and colons within the dis- 
tance a from the macroion surface are iV+ = 11 and 
N~ = 5, respectively. This means that the net charge 
of the entire complex, i.e. "macroion -f attached counte- 
rion -t- attached colons", is -|-44e. This is to be compared 
with the bare macroion charge of — 28e, which amounts 
to charge inversion of about 160% the original macroion 
charge. 




FIG. 1. The bird's-eye view of the screening atmosphere 
within 3a from the macroion under "standard" conditions. 
Macroion is the red ball in the middle. Multivalent coun- 
terions of valence Z = 7 and monovalent coions are shown 
in dark blue and light blue, respectively. Macroion radius 
-Ro = 3a, charge Qo = — 28e. Temperature is chosen such 
that Ta = 29.4. Note that significant condensation of coions 
on the counterions is observed. For this reason, correlations 
between adsorbed multivalent counterions are nowhere near 
ideal Wigner crystal while V ~ 137 (Eq.(j^) is very large. 

Figure ^ (a) shows the radial distributions of co- and 
counter-ions charges 

p.(r)=eZ./^fc;i^df7r, (8) 

where s means either co- or counter-ions, Zs is, accord- 
ingly, either —1 or Z] summation runs over all ions of 
the given sort s, Ygi is the position vector of ion i of the 
sort s, and fir is the solid angle of directions of vector r. 



These results are consistent with the conclusion of 160 % 
charge inversion. Indeed, the distribution of the counte- 
rions, denoted by open bars, is sharply peaked at r = Rq, 
while that of the colons (shaded bars) is broad and de- 
tached from the macroion surface. Although at this stage 
we do not formulate any rigorous algorithmic definition 
as to which counterions are close enough to the macroion 
to be called "bound," we note that the peak in the ra- 
dial density distribution of counterions is sharp enough 
to provide for quite clear distinction between bound and 
unbound ions. We therefore rely on this sharp peak, and 
in what follows we describe as bound those counterions 
which belong to this peak. 
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FIG. 2. Charge inversion under "standard conditions", as 
in Figure |l|. (a) The radial distribution function of the 
charge ps{r) (Eq.(^)) of counterions (open bars) and that of 
colons (shaded bars) as a function of the distance r from 
the macroion center, (b) The integrated charge distribu- 
tion Q(r) of counterions plus coions (Eq.(|^). The portion 
Q{r)/\Qo\ > 1 corresponds to charge inversion. 

Figure || (b) depicts the integrated charge of the mov- 
able ion species (counterions and coions) of Fig. ^ (a), 
starting at the surface of the macroion, 



dr' 



(9) 



The portion above the baseline Q/\Qo\ = 1 corresponds 
to the charge inversion (this applies to all the follow- 
ing figures). The net amount of inverted charge reaches 
160% for this run, as stated above, and the Q{r) pro- 
file relaxes to neutrality in a distance of approximately a 
few a, thus suggesting once again that a significant pop- 
ulation of coions reside on the outer sides of condensed 
counterions. Fluctuations of Q{r) for r ^ Rq reflect 



density fluctuations which are much amplified because 
of the volume factor 47rr^. On the other hand, we ob- 
serve a nearly neutral region Q/\Qo\ ~ 1 extending for 
the distance comparable to the Bjerrum length £b out- 
side the charge inversion layer. Few ions exist in this 
region. This shows establishment of enhanced order due 
to strong Coulomb interactions. 

The electrostatic potential drop across the charge dis- 
tribution peak corresponds to energy change eAip « 
1.2e^/ea, which is five times the thermal energy fc^T. 
This implies strong binding of counterions to the 
macroion and coions to the counterions. In other words, 
this manifests very strongly non-linear screening com- 
pared with Debye-Hiickel screening of weakly coupled 
cases. Of course, this is by no means surprising given 
the small value of A Eq.(|^), as mentioned above. 




FIG. 3. The bird's-eye view of the screening atmosphere 
within 3a from the "large" macroion: macroion radius 
Ro — 8a, charge Qo = — 28e. Temperature is adjusted such 
that e'^/eRnT = const is the same as in Fig.hl. This means 



that F Eq.(g) is the same here and in Fig.^, while Ta is greater 
here than in Flg|l| by a factor of 8/3. Accordingly, stronger 
binding of monovalent coions (dark blue) to multivalent coun- 
terions (light blue) is observed, and adsorbed counterions are 
less strongly correlated. 

Speaking about the dynamics of equilibration, it is in- 
teresting to note that the buildup of counterions on the 
macroion occurs fairly quickly, in about lOOw^"^, which 
is of the order of 100 picosec for the typical numerical 



values of parameters, as suggested in Section II B. This 
time is much shorter than overall relaxation time of the 
system, suggesting that equilibration of plasma further 
away from macroion occurs fairly slowly. It is appeal- 
ing to guess that this fast buildup of screening (and even 
over-charging) layer is connected with the fact of strongly 
non-linear correlated screening. 



2. Other regimes 

The charge inversion for the macroion with a large ra- 
dius i?o = 8a is depicted in Figs. Sand ^. Other param- 
eters are the same as those of Fig.^ except for the lower 
temperature (Fa = 78.4) to keep Fg = const Eq.(||). 
We again observe sparsely distributed counterions on the 
macroion surface. In this case, however, binding of the 
counterions to the macroion is loose, and their radial dis- 
tribution in Fig. ^ (a) is almost as broad as that of the 
colons. The counterion charge is better canceled on each 
site by the condensed colons than in Fig. ||. 

OiKH I ^ 




FIG. 4. Charge inversion of a "large macroion" - the sys- 
tem shown in Fig.|l]. Plot format and notations are the same 
as in Fig-^ Simulation parameters are also the same as those 
of Fig.hl, except that the temperature is adjusted to keep 
constant the quantity e^ /tRoT — const. Counterions are 
loosely bound to the macroion as Qo is fixed, and although 
the amount of inverted charge is smaller compared to that in 
Fig.^, it is still significant. 

We note that the number of condensed ions to the 
macroion surface in Fig. || is ^ 13 and ^ 66, 
where the number of is comparable to that in Fig. 
|l|. This is consistent with the fact that each coun- 
terion occupies, roughly, a neutralizing region on the 
macroion surface, similar to the Wigner-Seitz cell of 
Wigner crystal. With charge density of the macroion 
surface a = Qq/AttRq, the size of such neutralizing re- 
gion, or cell, is proportional to the size of the macroion: 
eZ = TTcrRl^, or R^s = 2Ra{Ze/\Qa\)'^/'^ . In other 
words, the neutralizing number of counterions {Rq/ R^s)^ 
stays unchanged as long as the macroion charge Qo is 
fixed. The inverted charge in Fig. ^ (b) is about 40% 
the original charge of the macroion, which is less than 
that in Fig. ||. The electrostatic potential drop across 



the macroion surface is consistently less than the ther- 
mal energy, eA^s 0.05e^/ea < ksT ^ 0.09e^/ea. The 
linear Debye-Hiickel theory nearly applies in this case. 

We found similar features, based on identification of 
bound ions in the peak of their radial distribution, also 
for the parameters further away from our standard con- 
ditions. For instance, we mention here in passing the 
case of the counterions with smaller valence Z = 3. For 
them, it takes somewhat less than 1 x lQ^uj~^ to reach the 
stationary state, and the attained peak height is lower, 
about 70% the macroion charge, as shown in Fig. ^. This 
will be discussed in greater details in one of the sections 
below. 
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FIG. 5. Dependence of inverted charge on the radius of 
macroion _Ro shown for the counterions with the valence Z = 3 
and 7. The charge of the macroion is Qo = — 28e, and the 
number of colons A*'" = 335 (or 336) corresponds to the den- 
sity n~ ~ 1 X 10~^a~'^. The ordinate is the maximum of 
the integrated charge Q{r) (Eq.(|)), i.e. Qp^ak = max((3(r)), 
normalized by the macroion charge |Qo|. Each data point is 
an average of three runs, and a vertical bar shows the range 
of time variations. 



B. Changing macroion properties and temperature 



In the following figures. Figs. |5[-|9|, the ordinate Qpeak 
is the maximum of the integrated charge of the counte- 
rions plus colons, Eq.(^. Each data point is an average 
of three runs, and a vertical bar shows the range of time 
variations and deviations among the runs. 

The dependence of charge inversion on the radius of the 
macroion is depicted in Fig. ^ For different values of the 
radius, temperature is adjusted accordingly to keep un- 
changed the value of Tq cx 1/{RqT) Eq.(||). The valence 
of the counterions is chosen either Z = 3 or 7. The num- 
ber of counterions is iV+ = 121 and = 52 for Z = 3 
and 7, respectively, which is large compared to \Qo\/Ze 
required for charge neutralization of the macroion. These 
parameters are chosen in such a way that the number of 
colons, which is determined by neutrality condition, is 
virtually fixed, being N~ = 335 for the Z = 3 case and 
= 336 for the Z = 7 case. This corresponds to 



Eq.(3) moderately changing between 0.3a ^ 0.8a, and A 
Eq.(7) changing between 0.02a ~ 2.8a. 

In Fig. 1^, the inverted charge reaches its maximum for 
the radius Ro « 3a irrespectively of the valence Z. It falls 
off rapidly both for smaller and larger radii, and becomes 
insensitive to the radius of the macroion for Ro/a 3> 1. 
The net amount of the inverted charge is about 70% of 
the bare macroion charge Qo for Z = 3; it increases up to 
150% of Qo for Z = 7. We find that the charge inversion 
reaches maximum also at virtually the same radius Ro w 
3a even for the smaller number of counterions = 15 
(Z = 7), or for larger macroion charge Qo = — 42e. 

It is not difficult to understand qualitatively why the 
charge inversion decreases at both small and large values 
of macroion radius Ro, reaching a maximum in between. 
When Ro gets very large, the lateral spacings between 
bound counterions become too long to maintain correla- 
tions between them; on the other hand, when Ro gets 
too small, the increased repulsion of the inverted charge 
from the macroion becomes dominant. 
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Z^e2 /eaT 

FIG. 6. Dependence of inverted charge on temperature 
shown for counterions with different valence Z, which is given 
by a master curve. The abscissa is the ratio of the Coulomb 
energy of the counterions to thermal energy, Ze^ /eaT. The 
radius of macroion is Ro — 3a, and the number of colons is 
kept nearly the same, A'^" ~ 335 or 336 for Z = 5 and 7. 



former due to larger Coulomb binding energy, and higher 
temperatures are better to suppress the latter due to en- 
hanced thermal motion. 

Figure]^ shows that charge inversion Qpeak/\Qo\ is in- 
sensitive to the charge content of the macroion Qo for 
fixed value of Tq = Qoe/eRoT (Fa = 4.2Z or 6Z). The 
number of counterions attached to the macroion surface 
is in the range 8 15 for |Qo| = (14 - 42)e and Z = 7, 
which is a few times that of the neutralizing number of 
counterions, \Qo\/Ze. 
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FIG. 7. Dependence of inverted charge on the macroion 
charge Qo- The radius of macroion is Ro — 3a, the valence 
of counterions is Z = 7. Temperature is adjusted to keep the 
binding energy Fq = \Qo\e/eRQT constant as Qo varies. 

We note in passing that the geometrical capacity of the 
surface, controlled by non-Coulomb short range forces 
is still very far from exhausted, AnRo/Tra'^ ^ 36. The 
regime of closed and almost closed packing of the bound 
spheres on the macroion is examined in the recent work 
pl| . Interestingly, the effective valence of the counteri- 
ons Zeff, which is the charge of the counterion minus 
that of the condensed ions, increases with the charge of 
the macroion; it is Z^ff ~ 0.25Z for Qo ~ — 14e and is 
Zeff ~ 0.4Z for Qo = -42e. 



The effect of temperature on charge inversion is shown 
m Fig. I In this fi gure, the abscissa is F^ = Ze^ /eaT 
Eq.(^. As the figure indicates, the inverted charge for 
different values of valence form a master curve when plot- 
ted against F^ . The charge inversion is maximized at the 
intermediate temperature corresponding to Fq 45, or 
Ze^/eRoT ~ 15 (i?o = 3a). The value of the Debye 
length is k, 0.6a for Z = 3 and k. a for Z = 7 . 
For the low temperature side. Fa ^ 100, the integrated 
charge distribution Q(r) is sharply peaked as that of Fig. 
H (b), while at the high temperature side, F^ ~ 10, this 
distribution Q(r) is rugged and fluctuates considerably 
with time. It is consistent with that the maximal charge 
inversion is achieved through competition of counterion 
attachment to the macroion and colon condensation on 
the counterions. Lower temperatures are favored for the 



C. Changing counterion properties 

The dependence of inverted charge on the valence of 
the counterions is depicted in Fig. ^ Here, the macroion 
charge and radius are Qo = — 28e and Ro = 3a, re- 
spectively, and temperature is fixed such that F^/Z = 
e^ /eaT = 4.2. It is emphasized that no charge inversion is 
observed for monovalent counterions. The amount of the 
inverted charge Qpeak increases with the valence, which 
is well scaled by Qpeak Z^/^ ^j. z < 5. The Z > 5 
part can be fit by Qpeak ^ Z. The inverted charge is also 
an increasing function of the number of counterions and 
colons, as seen by the difference of the two curves for two 
densities in the figure. 




^ N=10O 



FIG. 8. Amount of inverted charge Qpeak increasing mono- 
tonically with the valence of counterions Z. The inverted 
charge is well fit by Qpeak ~ Z^^'^ for ^ < 5, and Qpeak ~ Z 
for Z > 5. Note that charge inversion occurs only for multi- 
valent counterions, i.e. Z >2. 
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FIG. 9. Dependence of inverted charge on the ionic 
strength m — {Z^N^ + N~)/V, where V is the volume of 
the simulation cell. Charge neutrality of the system is main- 
tained. The guide curve is for ni < 0.02/a^, and ni 
for rii > 0.05/a"^. The macroion radius is Rq — 3a, charge 
Qo = — 28e, the valence of counterions Z = 7, and Fa = 29.4. 

The dependence of inverted charge on the ionic 
strength, n/ = {Z'^N+ + N-)/V, is shown in Fig. |, 
where V = 47r(i?|,^ — -f^o)/^ is the domain volume. The 
amount of inverted charge Qpeak /\Qo\ increases mono- 
tonically with the ionic strength. The functional form 
of the scaling changes at nj ^ O.OS/a"^, as shown by 
fitting curves. The ionic strength of a Ca,'^^ ion and 



neutralizing colons in every 10 A cube yields 0.048/a'^ 
for a = 2A. The scaling Qpeak ^ n-j for the low 
ionic strength nj < 0.02/0*^ smoothly joins a linear scal- 
ing Qpeak ni for high ionic strength nj > 0.05/a'^. 
The non-dimensional parameter of the theory C = 
{R^sh-sf = l2aTaNei{e/\Qo\){Rl/Rlj), is calculated to 
be 0.7 for n/ ^ O.Ol/a"^ and Z = 7. The theory expects 
qm _ {NeiZf^ for C < 1, and Q^*'') ~ iV„Z for 
C ^ 1. The present simulation results agree with this 
theoretical prediction. 
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FIG. 10. Potential energy shown as a function of Ze^ /eaT 
for the runs with Z = 1 (cf. Fig. ^). The filled and open 
circles correspond to the potential energy of interaction of a 
macroion with counterions and colons, respectively, the tri- 
angles those between counterions and colons, and the squares 
the total potential energy. 



D. Measuring potential energy 

The potential energy presented in Fig. ^ is in line 
with the tendency of charge inversion dependence on vari- 
ations of the Coulomb coupling parameter (cf. Fig. |^). 
The potential energy for the interactions between coun- 
terions and the macroion (solid circles) is negative (at- 
tractive) and is minimized at the intermediate value of 
Fq — Ze^ /eaT ~ 50 where largest charge inversion is ob- 
tained. The potential energy of interactions between the 
counterions and colons, depicted by triangles, decreases 
remarkably with the increase in Fq. This corresponds 
to massive condensation of colons onto the counterions 
(similar to Manning-Onsager condensation) at low tem- 
peratures. This reduces the effective valence of the coun- 
terions, and the binding of counterions on the macroion 
surface is weakened, which tends to suppress the charge 
inversion. Thus, charge inversion becomes largest at the 
intermediate value of Ze^/eaT, as stated above. The 
potential energy of interaction between colons and the 
macroion (open circles) is positive (repulsive), and is 
maximized where the colons are closely located with the 
macroion by condensation to the counterions, namely at 
Fq ~ 50. On the other hand, the total potential energy 



(squares) decreases with the increase in the couphng pa- 
rameter Ta- 



current study. We also thank him for reading the paper 
and giving them valuable comments. 



IV. CONCLUSION 



In this paper, we showed the occurrence of giant charge 
inversion and examined its parameter dependences with 
the use of molecular dynamics simulations. The charge 
inversion was found to be based on the strong correlations 
of the multivalent counterions and colons, particularly on 
the surface of the macroion. Specifically, charge inversion 
was observed under the conditions for which the Coulomb 
coupling parameter was significantly larger than unity, 
r ^ 1. At the same time, charge inversion occurred 
only in the presence of multivalent counterions with Z > 
2. The counterions were attached to the surface of the 
macroion, while monovalent colons tended to condense on 
the counterions. This condensation, similar to Bjerrum 
pairing, is therefore identified as the factor limiting the 
amount of charge inversion. The amount of the inverted 
charge Qpeak was maximal at rather small radius of the 
macroion, and leveled off when radius becomes larger. It 
scaled linearly with the charge of the macroion Qq, and 
the ratio Qpeak/\Qo\ was independent of the macroion 
charge. 

With respect to the valence Z and the ionic strength 
nj = {Z'^N^ + N^)/V, the amount of inverted charge 
scaled as Q ~ (ZniY^^ for the valence Z < 5 or n/ < 
0.02/a^. As noted in Sec. Ill C, this ionic strength corre- 
sponds to a Ca^''" ion in every lOA cube. The inverted 
charge scaled as Q Znj for Z > 5 or n/ > 0.05/a'^. 
This agreed with the theory of giant charge inversion 
The net inverted charge of nearly up to 160% the 
bare charge of the macroion was achieved at the medium 
temperature Ze^ /eRoT ^ 15, due to the competition of 
multivalent counterion attachment to the macroion and 
monovalent colon condensation on the counterions; the 
former was stronger at lower temperatures, and the lat- 
ter was suppressed at higher temperatures. 

In the present study, the macroion was assumed to be 
immobile. From the application points of view, it might 
be informative to study the distribution of counterions 
and colons around a moving macroion and also the effect 
of an applied electric field. The study of such cases is 
reported in a separate paper. The results indicate that 
a formed complex of a macroion and counterions drifts 
along the electric field in the direction implied by the 
inverted charge, and that charge inversion is not altered 
until the electric field exceeds a critical value |27| . 
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